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ABSTRACT 



We analyse the performance of twelve different implementations of Smoothed Par- 
ticle Hydrodynamics (SPH) using seven tests designed to isolate key hydrodynamic 
elements of cosmological simulations which are known to cause the SPH algorithm 
problems. In order, we consider a shock tube, spherical adiabatic collapse, cooling 
flow model, drag, a cosmological simulation, rotating cloud-collapse and disc stability. 
In the implementations special attention is given to the way in which force symmetry 
is enforced in the equations of motion. We study in detail how the hydrodynamics 
are affected by different implementations of the artificial viscosity including those 
with a shear-correction modification. We present an improved first-order smoothing- 
length update algorithm that is designed to remove instabilities that are present in 
the Hernquist and Katz ( 1989| ) algorithm. Gravity is calculated using the Adaptive 



Particle-Particle, Particle-Mesh algorithm. 

For all tests we find that the artificial viscosity is the single most important 
factor distinguishing the results from the various implementations. The shock tube and 
adiabatic collapse problems show that the artificial viscosity used in the current HYDRA 
code performs relatively poorly for simulations involving strong shocks when compared 
to a more standard artificial viscosity. The shear-correction term is shown to reduce the 
shock-capturing ability of the algorithm and to lead to a spurious increase in angular 
momentum in the rotating cloud-collapse problem. For the disc stability test, the 
shear-corrected and current hydra artificial viscosities are shown to reduce outward 
angular momentum transport. The cosmological simulations produce comparatively 
similar results, with the fraction of gas in the hot and cold phases varying by less than 
10% amongst the versions. Similarly, the drag test shows little systematic variation 
amongst versions. The cooling flow tests sho w tha t implementations using the force 



symmetrization of Thomas and Couchman ( 1992 ) are more prone to accelerate the 
overcooling instability of SPH, although the problem is generic to SPH. The second 
most important factor in code performance is the way force symmetry is achieved in 
the equation of motion. Most results favour a kernel symmetrization approach. The 
exact method by which SPH pressure forces are included in the equation of motion 
appears to have comparatively little effect on th e resu lts. Combining the equation 
of motion presented in Thom as and Couchman ( 199S ) with a modification of the 
Monaghan and Gingold ( 1983| ) artificial viscosity leads to an SPH scheme that is both 



fast and reliable. 
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Smoothed Particle Hydr odynamics (SPH) (Gingold & Mon- 



aghan 1977; Lucy 1977) is a popular numerical technique 
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for solving gas-dynamical equations. SPH is unique among 
numerical methods in that many algebraically equivalent - 
but formally different - equations of motion may be derived. 
In this paper we report the results of a comparison of sev- 
eral implementations of SPH in tests which model physical 
scenarios that occur in hierarchical clustering cosmology. 

SPH is fundamentally Lagrangian and fits well with 
gravity solvers that use tree structures (Hernquist & Katz 
1989, hereafter HK89) and mesh methods supplemented 
by short range forces (Evrard 1988; Couchman, Tho mas fc 
Pea rce 1 995, hereafter CTP95). In an adaptive form (Wood 
198l|)7the algorithm lends itself readily to the wide range 
of densities encountered in cosmology, contrary to Eulerian 
methods which require the storage and evaluation of numer- 
ous sub-grids to achieve a similar dynamic range. SPH also 
exhibits less numerical diffusivity than comparable Eulerian 
techniques, and is much easier to implement in three dimen- 
sions, typically requiring 1000 or fewer lines of FORTRAN 
code. 

The main drawback of SPH is its limited ability to fol- 
low steep density gradients and to correctly model shocks. 
Shock- capturing requires the introdu ction of an artificial vis- 
cosity (Monaghan & Gingold 1983). A number of different 



putational efficiency. This is a guiding principle in our in- 
vestigation. 

The seven tests used in this study are: 



alternatives may be chosen and it is not clear whether one 
method is to be preferred over another. The presence of 
shear in the flow further complicates this question. 

Much emphasis has been placed upon the performance 
of SPH with a s mall number of particles (order 100 or fewer). 
Initial studies ( Evrard 1988[ ) of SPH on spherical cloud col- 
lapse i ndicated acceptable performance when compared to 



low resolution Eulerian simulations, with global properties, 
such as total thermal energy, being reproduced well. A more 
recent study (Steinmetz & Muller 1993, hereafter SM93), 
which compares SPH to mod ern Eulerian techniques (th e 
Piecewise Parabolic Method (Collela & Woodward 1984), 



and Flux-Corrected- Transport methods) has shown that the 
performance of SPH is not as good as initially believed, 
and that accurate reproduction of local physical phenom- 
ena, such as the velocity field, requires as many as 10 4 parti- 
cles. In the context of cosmology with hierarchical structure 
formation, the small- N performance remains critical as the 
first objects to form consist of tens - hundreds at most - of 
particles and form, by definition, at the limit of resolution. 
It is therefore of crucial importance to ascertain the per- 
formance of different SPH implementations in the small- N 
regime. Awareness of this has caused a number of authors 
to perf orm deta i led tests on the limits of SPH ( Dwens & 
Villum jen~1997| ; [Bate fc Burkert 199t| ). To address these 
concerns, some of the tests we present here are specifically 
designed to highlight differences in performance for small N. 

Our goal in this paper is to detail systematic trends in 
the results for different SPH implementations. Since these 
are most likely to be visible in the small- regime we con- 
centrate on smaller simulations, using a larger number of 
particles to test for convergence. Because of the importance 
of the adaptive smoothing length in determining the local 
resolution we pay particular attention to the way in which it 
is calculated and updated. We are also concerned with the 
efficiency of the algorithm. Realistic hydrodynamic simula- 
tions of cosmological structure formation typically require 
10 4 or more time-steps. Therefore when choosing an imple- 
mentation one must carefully weigh accuracy against com- 



Sod shock (section 3.1) 



Spherical collapse (section 3.2) 



Cooling near density jumps (section 3.3) 
Drag on a cold clump(section 3.4) 



• Hierarchical structure formation (section 3.5) 

• Disc formation (section 

• Disc stability (section 

Each of these tests is described in the indicated sec- 
tion. The tests investigate various aspects of the SPH al- 
gorithm ranging from explicit tests of the hydrodynamics 
to investigations specific to cosmological contexts. The Sod 
shock (1978), although a relatively simple shock configu- 
ration, represents the minimum flow discontinuity that a 
hydrodynamic co de should be able to reproduce. The spher- 
ical collapse test (Evrard 198S ) , although idealised, permits 
an assessment of the resolution necessary to approximate 
spherical collapse. It also allows a comparison with other 
authors' results and with a high resolution spherically sym- 
metric solution. The remaining tests are more closely tied 
to the arena of cosmological simulations. The cooling test 
looks at the problems associated with modelling different 
gas phases with SPH. A cold dense knot of particles embed- 
ded in a hot halo will tend to promote cooling of the hot gas 
because of the inability of SPH to separate the phases. The 
drag test looks at the behaviour of infalling satel lites and 
the over -merging problem seen in SPH simulations ( Frenk ei 



/. 1996). Finally, three tests consider the ability of the SPH 



algorithms to successfully model cosmic structure. First we 
look at the overall distribution of hot and cold gas in a 
hierarchical cosmological simulation, followed by an inves- 
tigation of disc formation from the collapse of a rotating 
cloud (Navarro & White 1993), and the transport of angular 
momentum in discs. In each case we compare the different 
algorithms and assess the reliability of the SPH method in 
performing that aspect of cosmological structure formation. 

The layout of the paper is as follows. Section ^| reviews 
the basic SPH framework that we use, including a descrip- 
tion of the new approach developed to update the smoothing 
length. Next we examine the equations of motion and inter- 
nal energy, and discuss the procedure for symmetrization of 
particle forces. In section ^ we present the test cases as listed 
above. Each of the subsections is self contained and contains 
a description and motivation for the test together with re- 
sults and a summary comparing the relative merits of the 
different implementations together with an assessment of the 
success with which the SPH method can perform the test. 
Section ^ briefly summarises the main overall conclusions 
to be drawn from the test suite, indicating where each im- 
plementation has strengths and weaknesses and makes rec- 
ommendations for the implementation which may be most 
useful in cosmological investigations. 

During final prepartion of this paper, a preprint detail- 
ing a sim ilar investigation was circulated by Lombardi et 
al. jl998|). 



2 IMPLEMENTATIONS OF SPH 
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2.1 Features common to all implementations 

All of the implementations use an adapti ve particle-particle - 
particle-mesh (AP 3 M) gravity solver (Couchman 1991). 
AP 3 M is more efficient than standard P 3 M, as high density 
regions, where the particle-particle summation dominates 
calculation time, are evaluated using a further P 3 M cycle 
calculated on a high resolution mesh. We denote the process 
of placing a high resolution mesh 'refinement placing' and 
the sub- meshes are termed 'refinements'. The algorithm is 
highly efficient with the cycle time typically slowing by a fac- 
tor of three under clustering. The most significant drawback 
of AP 3 M is that it does not yet allow multiple time-steps. 
However, the calculation speed of the global solution, com- 
pared with alternative methods such as the tree-code, more 
than outweighs this deficiency. Full details of the adaptive 
scheme, in particular ac curacy and timing information, may 
be found in Couchman (|l99l|) and CTP95. 

Time-stepping is performed using a Predict-Evaluate- 
Correct (PEC) scheme. This scheme is tested in detail 
against leapfrog and Runge-Kutta methods in CTP95. The 
value of the time-step, dt, is found by searching the par- 
ticle lists to establish the Courant conditions of the ac- 
celeration, dt a , and velocity arrays, dt v . In this paper we 
introduce a further time-step criterion, dth, which pre- 
vents particles travelling too far within their smoothing 
radius, (see section |2,2| ). dt is then calculated from dt < 
k mm(0Adt v ,0.25dt a ,0.2dth) where s is a normalisation 
constant that is taken equal to unity in adiabatic simula- 
tions. In simulations with cooling, large density contrasts 
can develop, and a smaller value of k is sometimes required. 
We do not have a Courant condition for cooling since we im- 
plement it assuming constant density (see below and CTP95 
for further details). 

SPH uses a 'smoothing' kernel to interpolate local hy- 
drodynamic quantities from a sample of neighbouring points 
(particles). For a continuous system an estimate of a hydro- 
dynamic scalar A(r) is given by, 



(A(r))= d A r'A{r')W{r -r',h), 



(1) 



where h is the 'smoothing length' which sets the maximum 
smoothing radius and W(r,h), the smoothing kernel, is a 
function of |r|. For a finite number of neighbour particles 
the approximation to this is, 



(Mr)) 



Mjj) 
p( r j) 



W(r - 



,h), 



(2) 



where the radius of the smoothing kernel is set by 2h (for 
a kernel with compact sup port). The smoothing kernel u sed 
is the so-called i?2-spline (Monaghan & Lattanzio 1985), 



W(r,h) = 



Ws(r/h) 



h 3 

where if x = r/h, 



1 



4-6x 2 +3x 3 , < x < 1; 



1 < x < 2; 
x > 2. 



(3) 



(4) 



The kernel gradient is modified to give a small re pulsive 
force for close particles ( Thomas fc Couchman 1992 ), 



dWs(x) 
dx 




0<x< 2/3; 
3x), 2/3<x< 1; 
c) 2 , 1< x < 2; 
x > 2. 



(•>) 



The primary reason for having a non-zero gradient at the 
origin is to avoid the artificial clustering noted by Schiissler 
& Schmidt (1981). (Some secondary benefits are discussed 
by Steinmetz 1996.) 

In standard SPH the value of the smoothing length, 
h, is a constant for all particles resulting in fixed spatial 
resolution. Fixed h also leads to a slow-down in the calcu- 
lation time when particles become clustered - successively 
more particles on average fall within a particle's smoothing 
length. In the adaptive form of SPH (Wood 1981) the value 



of h is varied so that all particles have a constant (or ap- 
proximately constant) number of neighbours. This leads to 
a resolution scale dependent upon the local number density 
of particles. It also removes the slow-down in calculation 
time since the number of neighbours is held constant pro- 
vided that the near neighbours can be found efficiently. In 
this work we attempt to smooth over 52 neighbour particles, 
and tests on the update algorithm we use (see section 2.z) 
show that this leads to a particle having between 30 and 80 
neighbours, whilst the average remains close to 50. 

A minimum value of h is set by requiring that the SPH 
resolution not fall below that of the gravity solver. We de- 
fine the gravitational resolution to b e twice the S2 softening 
length (Hockney & Eastwood 1981), e, as at this radius the 
force is closely equivalent to the 1/r 2 law. We quote the 
equivalent Plummer softening throughout the paper, as this 
is the most common force softening shape. Since the min- 
imum resolution of the SPH kernel is the diameter of the 
smallest smoothing sphere, 4h m in , equating this to the min- 
imum gravitational resolution yields, 



2e = AK 



(6) 



Unlike other authors, once this minimum h is reached by 
a particle, smoothing occurs over all neighbouring particles 
within a radius of 2h m i n . As a result setting a lower h m in 
increases efficiency as fewer particles contribute to the sam- 
pling, but this leads to a mismatch between hydrodynamic 
and gravitational resolution scales which we argue is unde- 
sirable. 

When required, radiative cooling is implemented in an 
i ntegral form that assumes con stant density over a time-step 
(Thomas & Couchman 1992). The change in the specific 



energy, e, is evaluated from. 



de 

X 



-At, 



Pi 



(7) 



where At is the time-step, m the number density and n 2 A 
is the power radiated per unit volume. In doing this we cir- 
cumvent having a Courant condition for cooling, and hence 
it never limits the time-step. 



2.2 An improved first-order smoothing-length 
update algorithm 

As stated earlier, in the adaptive implementation of SPH the 
smoothing length, h, is updated each time-step so that the 
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Figure 1. Improve men t in the neighbour counts as each component of the new algorithm is added. A rotating cloud collapse 
problem (see section |3,6| ) was repeated with each of our modifications to the HK89 algorithm. The improvement from the first panel 
(upper left) to the third panel (bottom left) is clear, the final panel shows a slight reduction in the range of neighbours and in the 
oscillation of the counts. The large increase in the maximum number of neighbours is due to the h m i n limit being reached. Time 
units are given in Gyr. 



number of neighbours is held close to, or exactly at, a con- 
stant. Ensuring an exactly constant number of neighbours is 
computationally expensive (requiring additional neighbour- 
list searching) and hence many researchers prefer to update 
h using an algorithm that is closely linked to the local den- 
sity of a particle. 

We have used two guiding principles in the design of 
our new algorithm. First, since gravitational forces are at- 
tractive the algorithm will spend most of its time decreasing 
h (void evolution is 'undynamic' and poses little problem to 
all algorithms). Second, smoothing over too many particles 
is generally preferable to smoothing over too few. Despite 
some loss of spatial resolution and computational efficiency, 
it does not 'break' the SPH algorithm. Smoothing over too 
few particles can lead to unphysical shocking. 

A popular method for updating the smoothing length 
is to predict h at the next time-step from an average of the 
current h and the h implied by the number of neighbours 
found at the current time-step (HK89). This is expressed as, 



hi 



1 + 



N 3 



NT 



1/3 



(8) 



where h™ is the smoothing length at time-step n for particle 
i, N s is the desired number of neighbours and N™ -1 is the 
number of neighbours found at step n — 1. The performance 
of this algo rithm for a rotating cloud collapse problem (see 
section 3i3) is shown in the top left panel of Fig. [j]. The ro- 
tating cloud collapse is a difficult problem for the h update 
algorithm to follow since the geometry of the cloud rapidly 
changes from three to two dimensions. The plot shows that 
both the maximum and minimum number of neighbours (se- 
lected from all the SPH particles in the simulation) exhibit 
a significant amount of scatter. Clearly, the maximum num- 
ber of neighbours increases rapidly once the h = h m i n limit 
is reached. 

The rotating cloud problem demonstrates how this al- 
gorithm may become unstable when a particle approaches a 
high density region (the algorithm is quite stable where the 
density gradients are small). If the current h is too large, 
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the neighbour count will be too high leading to an under- 
estimate of h. At the next step too few neighbours may be 
found. The result is an oscillation in the estimates of h and 
in the number of neighbours as a particle accretes on to 
a high density region from a low density one. This oscilla- 
tory behaviour is visible in the fluctuations of the minimum 
number of neighbours in top left panel of Fig. |lj 

The instability can be partially cured by using an ex- 
tremely small time-step. However, this is not practical as 
simulations currently take t housan ds of time-steps. A solu- 
tion suggested by Wadsley ( 1997 ) that helps alleviate the 
discontinuity in the number of neighbours, is to 'count' 
neighbours with a smoothed kernel rather than the usual 
tophat. This (normalised) weighted neighbour count is then 
used in the h update equation rather than N™~ . The insta- 
bility in the standard HK89 algorithm is caused by the sharp 
discontinuity in the tophat at r — 2h. Hence we require that 
the new kernel smoothly approach at 2h. Secondly, we 
wish to count the most local particles within the smooth- 
ing radius at full weight. Hence we consider a kernel that is 
unity to a certain radius followed by a smooth monotonic 
decrease to zero. We choose, 

W (rlh)-i l > 0<r/h<3/2- 

Wnn[ / ) -\nW s (Mr/h-3/2)), 3/2 < r/h < 2, w 

where W s (x) is the normalised _B2-spline kernel. We have ex- 
perimented with changing the radius at which the counting 
kernel switches over to the spline, and a value of r = 3h/2 
has proven to be optimal, providing a good balance between 
the smoothness of the variation of the estimate and the close- 
ness of actual number of neighbours to the desired number. 
At this value approximately half of the smoothing volume 
is counted at full weight. For smaller values the smoothed 
estimate becomes progressively more unreliable. Conversely 
as the limiting radius is increased the gradient of the kernel 
at r ~ 2h becomes too steep. The improvement in fluctua- 
tion of the maximum and minimum number of neighbours 
can be seen in the top right panel of Fig. [j]. 

The next step in the construction of the new algorithm 
is to adjust the the average used to update h. The pri- 
mary reason for doing this is to avoid sudden changes in 
the smoothing length. Whilst the neighbour counting kernel 
helps to alleviate this problem, it does not remove it entirely. 
Secondly, since we shall limit the time-step by only allowing 
the particles to move a certain fraction of h it is also useful 
to limit the change in h. The motivation here is that a parti- 
cle which is approaching a high density region, for example, 
and is restricted to move 0.2h per time-step should only be 
allowed to have h — > 0.9h (since the smoothing radius is 2h). 
However, a particle at the centre of homologous flow must 
be able to update faster since collapse occurs from all direc- 
tions. To account for this it is helpful to permit a slightly 
larger change in h, h — > 0.8h for example. 

Setting s = (Ns/N™^ 1 ) 1 / 3 , we may express equation ^ 

as 



h™ — h™ 1 (1 — a + as), 



(10) 



where a is a weighting coefficient, and for equation g| 
a = 0.5. We have tested the performance of this average for 
a 6 [0.2,0.5]. A value of a = 0.4 proved optimal, reducing 
scatter significantly yet allowing a sufficiently large change 
in h. A problem remains that if s ~ then h™ = 0.6/i™~ , 
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0.8+0.2s weighted average 
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Figure 2. Weighting function compared to simple weighted av- 



which represents a large change if one limits the time-step 
according to an h/v criterion. Hence we have implemented 
a scheme that has the asymptotic property h™ = 0.8/i™~ , 
but for small changes in h it yields h" — h%~ 1 (0.6 + 0.4s). 
The function we use for determining the weighting variable 
a is, 



0.2(1 + s 2 ), 
0.2(1 + l/s 3 ), 



s < 1; 
s > 1. 



(11) 



A plot of this function compared to the 0.6,0.4 weighted 
average can be seen in Fig. ^. The lower left panel of Fig. [j] 
shows that introduction of this average reduces the scatter 
in both the maximum and minimum number of neighbours. 

Even with the improvements made so far it remains 
possible that the particle may move very quickly on to a 
high density region causing a sudden change that cannot be 
captured by the neighbour counting kernel. Thus it is sensi- 
ble to limit the time-step according to a 'Courant condition' 
which is the smoothing radius divided by the highest veloc- 
ity within the smoothing radius. Further, the velocity v r , 
must be measured within the frame of the particle under 
consideration, v r = Vi — Vj. The new time-step criterion 
may be summarised, 

dt — C'h min(/ii/ ma,x(\vi — Vj |), (12) 

» i 

where Ch is the Courant number and the i, j subscripts de- 
note a reduction over the variable indicated. We have chosen 
Ch = 0.2, which is usually the limiting factor in time-step se- 
lection. The gains in introducing this condition are marginal 
for the neighbour count in this test (maximum and minimum 
values are within a tighter range and show slightly less os- 
cillation) . 

When taken together all of these adjustments combine 
to make a scheme that is both fast and very stable. The final 
result of including these improvements is seen in the lower 
right hand panel of Fig. |l|. 

2.3 Equations of motion 

The SPH equation of motion is derived from, 
dv 1, 



dt 



-VP, 



P 



(13) 



using identities involving the pressure and density. An ex- 
cellent review of this derivation, and why so many different 
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schemes are possible, may be found in Monaghan (1992). 
In short, different identities produce different equations of 
motion and a clear demonstration of this can be seen by com- 
paring the equations of motion of HK89 to those of SM93. 

Once an adaptive scheme is implemented in SPH, neigh- 
bour smoothing develops a dualism - the smoothing may be 
interpreted as either a "gather" or a "scatter" process. If 
we attempt to smooth a field at the position of particle i, 
then the contribution of particle j to the smoothed estimate 
may be evaluated using the value of h from either particle 
i (gather) or particle j (scatter). For constant h SPH the 
schemes are equivalent. 

The uncertainty in the neighbour smoothing may be 
circumvented by using a form alism based upon the av erage 
of the two smoothing lengths ( Evrard 1988 ; Benz 199C ). Un- 
der this prescription, the gather and scatter interpretations 
are equivalent since the smoothing length for particles i and 
j are the same. We call this prescription /i-averaging. The 
density estimate under this prescription is given by, 



(P(n)) 



where, 



JV 

£ 

j=i 



mjW(ri 



hij = (hi + hj)/2, 



(14) 



(15) 



and the neighbour search is conducted over particles for 
which m < 2hij. 

Most h- averaging schemes use the arithmetic mean of 
the smoothing lengths but it is possible to consider other av- 
erages (any 'average' that is symmetric in i — j is potentially 
acceptable). Two other averages of interest are the geomet- 
ric mean and the harmonic mean. Note that the harmonic 
and geometric means both pass through the origin where 
as the arithmetic mean does not. This has potentially im- 
portant consequences when particles with large h interact 
with particles that have a small h. This situation occurs at 
the boundary of high density regions in simulations with ra- 
diative cooling (see section 3.3 for a detailed discussion of 
this). 

An alternative way of circumventing the smoothing un- 
certainty (HK89) is to combine the "gather" and "scatter" 
interpretations into one hybrid framework by averaging the 
kernels. In this scheme the density estimate is given by, 

N 

(P(n)) =^2m j [W(r i -r j ,h i ) + W(r i -r j ,h j ))/2, (16) 

and we denote this scheme kernel averaging. The averaged 
kernel then replaces the normal kernel in all equations. 
Rewriting equation Il3l using the identity, 



VP P P„ 

— = v- + — Vp, 

p p p 



(17) 



the SPH equation of motion with kernel averaging becomes, 



dvi 
~dt 



Pi , Pi 



VilWin-r^hJ + Wtn-r^hj)]^. 



(18) 



This equation of motion is used in SM93. 



Thomas & Couchman (1992, hereafter TC92), present 
an alternative prescription where the force term is sym- 
metrized and the density is calculated under the gather in- 
terpretation. Using the standard pressure and density iden- 
tity in equation 17], the acceleration for particle i is written, 



dvi 



fi 



N 

mi 

3=1 



N 

£< 

3=1 



P 

1 2 
Pi 



ViW(r t ~r J ,h l ) 



N 



3 = 1 



(19) 



This symmetrization is used in the current implementation 
of the publically available code, hydra. When combined 
with an artificial viscosity that does not require the pre- 
computation of all density values this scheme is extremely 
efficient. 

In deriving equation ^ the approximation, 



ViW(n -rj,hi) 



(20) 



was used. This approach to symmetrization of the equa- 
tion of motion is fundamentally different to the other two 
schemes which both result in a kernel symmetric in i and 
j. The approximation is not correct to first order in h and 
may introduce small errors. If this symmetrization is sup- 
plemented by either /i-averaging or kernel averaging (there 
is no argument against this) then the substitution is correct 
to first order in h (but not in V/i). 

When SPH is implemented in its adaptive form equa- 
tion |l8| will not guarantee conservation of momentum. To do 
so one must alter the neighbour list so that for any particle i 
which smooths over particle j the converse is also true. This 
will not necessarily be true if one smooths over the nearest 
52 neighbours for each particle. Note that equation |l^ does 
conserve momentum regardless of the neighbour list used. 

In the SPH-AP 3 M solver an SPH particle only has 
gas forces evaluated within a refinement if it has a suffi- 
ciently small h (see CTP95 for a detailed discussion of this 
topic). The introduction of /i-averaging complicates refine- 
ment loading since a particle no longer has a well defined h 
(one must consider all hij values). To establish which parti- 
cles are to be evaluated, one must calculate all the hij values 
for a particle, which is computationally costly. As a com- 
promise, we place particles in a refinement only if 1.4 x 2h 
is smaller than the sub-mesh search length (note that for 
arithmetic ft-averaging a search to 2 X 2h will ensure that 
all neig hbours a re fou nd and we use this search radius in 
sections 3~5| , 3.6 & 3.7). This does not absolutely guarantee 
that all the neighbour data for a particle will be placed in 
the refinement data arrays but we have found no detectable 
difference between this method and the full hij method. 



2.4 Internal energy equation 

The SPH internal energy equation is derived from, 

du P„ ,„„. 

— = V.w. (21) 

dt p v ' 

The SPH estimate for V .v may be used to calculate du/dt 
directly, yielding the internal energy equation for particle i, 

dei Pi , 



dt 



-V.Vi 



(22) 

© 1998 RAS, MNRAS 000, 



Comparative study of SPH implementations 7 



Explicitly writing the summation for the SPH estimate 
gives, 

d£j \ - Pi 

dt 



/ Jo 

V Pi 



,h). 



(23) 



Substituting h = hij and inser ting the artificial viscosity 
Pi — *• Pi + pfllij/2 (see equation |2q ), yields the internal en- 
ergy equation used in SM93. Whilst not strictly compatible 
with equation [iS] (equation ^ has no dependence upon Pj) 
i t can never theless be shown that energy will be conserved 
( |Benz 199CD - 

An internal energy equation may be constructed using 
the same symmetrization as equation [Hj, yielding, 



dti 
~dt 



1 



■Vj).-^ViW(ri 



3=1 

JV 
3=1 



Pi 



,hi) 



rrij (vi 



(24) 



A similar equation was considered in TC92, and was shown 
to exhibit excellent energy conservation. In CTP95 the en- 
tropy scatter produced by this equation was compared to 
that of equation |22| and was shown to be significantly larger. 
Therefore to avoid this problem we only consider the energy 
equation 

It has been shown by N elson & Papaloizou (1993,1994) 
and Serna, Alimi & Chieze (1996), that inclusion of the V/i 
terms, in both the equation of motion and internal energy, is 
necessary to ensure conservation of both energy and entropy. 
We chose to neglect these terms due to the overhead involved 
in computing them. 



2.5 Artificial viscosity algorithms 

An artificial viscosity is necessary in SPH to dissipate con- 
vergent motion and hence prevent interpenetration of gas 
clouds (|Monaghan fc Gingold 1983). 



We investigate four different types of artificial viscosity. 
The first considered is the implementation of TC92, where 
an additional component is added to a particle's pressure, 



Pi + pi[~aahiV.Vi + l3(hiV.Vi) 2 }, V.Wi < 0; 



Pi, 



V.m > 0, 



(25) 



where a is the sound speed for the particle. Typically, the 
viscosity coefficients are a — 1 and (3 = 2. This algorithm 
is a combination of the Von Neumann-Richtmeyer and bulk 
viscosity prescriptions (see Gingold & Monaghan 1983). A 
notable feature is that it uses a V.Vi 'trigger' so that it only 
applies to particles for which the local velocity field is con- 
vergent. This is different from most other implementations, 
which use an rij.Vij < trigger. With this artificial viscos- 
ity the first term in equation [l9] does not depend upon the 
density of particle j. This is advantageous numerically since 
it is not necessary to construct the neighbour lists twice. In 
this circumstance one can reduce the SPH search over parti- 
cles to a primary loop, during which the density is calculated 
and the neighbour list is formed and stored, followed by a 
secondary loop through the stored neighbour list to calculate 
the forces and internal energy. 

A modification of this artificial viscosity prescription is 



to estimate the velocity divergence over a smaller number of 
neighbours found by searching to a reduced radius. This was 
motivated by the observation that in some circumstances the 
gas does not shock as effectively as when a pairwise viscosity 
is employed. The reduced search radius leads to a higher 
resolution (but likely noisier) estimate of V.v. 

The third artificial viscosit y we consider is the standar d 
'Monaghan' artificial viscosity (Monaghan & Gingold 1983). 



This artificial viscosity has an explicit i — j particle label 
symmetry which is motivated to fit with equation [l8| In 
this algorithm the summation of P/p 2 terms in equation |li| 
is extended to include a term fly , which is given by, 



n, 



where, 



-ctfiijdj + (3p\ 

Pi] 



hijVij.rij/ir^ + u 2 ), 
0, 



<0; 
>0, 



(26) 



(27) 



where the bar denotes arithmetic averaging of the quantities 
for particles i and j and v 2 = O.Ol/i^ is a term included to 
prevent numerical divergences. Again, typical values for the 
coefficients are a = 1 and (5 = 2 although for problems in- 
volving weak shocks values of 0.5 and 1, respectively, may be 
preferable. This artificial viscosity has the same functional 
dependence upon h as the V.u version. Since this algorithm 
utilises a pairwise evaluation of the relative convergence of 
particles it is capable of preventing interpenetration very ef- 
fectively. A drawback is that it requires the neighbour list 
for particles to be calculated twice since the force on particle 
i depends explicitly on the density of particle j (storing all 
the neighbour lists is possible, but in practice too memory 
consuming) . 

One major concern about this algorithm relates to its 
damping effect on angular momentum. In the presence of 
shear flows, V.i> = 0, V x v > the pairwise rij.Vij term 
can be non zero and hence the viscosity does not vanish. 
This leads to a large shear viscosity which is highly unde- 
sirable in simulations of disc formation, for example. A way 
around this problem i s to add a sh ear-correcting term to 
the artificial viscosity (Balsara 1995). The modification we 



consider is given by Steinmetz (1996); 



Hij 
where, 

Hi — 
and, 

/* = 



Hij fi; 



fi + f] 



(28) 



(29) 



(30) 



(V.»>. | + | (V x v) i | + O.OOOlft/h, ' 

For pure compressional flows / = 1 and the contribution of 
Hij is unaffected whilst in shear flows f — and the viscos- 
ity is zero. This m odification has been studied by Navarro & 
Steinmetz (1997) who found that the dissipation of angular 



momentum is drastically reduced in small- N problems using 
this method. However, of concern is whether this modifica- 
tion leads to poorer shock-capturing. This may arise due 
to sampling error in the SPH estimates of the velocity di- 
vergence and curl. We compare schemes using the shear- 
corrected viscosity against the standard Monaghan viscos- 
ity. 
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The efficiency gained from having an artificial viscosity 
that does not depend upon the density of particle j has 
motivated us to consider the following modification of the 
Monaghan viscosity: we use the same quantities as equation 
Eq, except that, 



Pii -» Pij = + (hi/hj) 3 )/2, 



(31) 



which provides an estimate of pij. The estimate is based 
upon the approximation pj ~ pi /if /ft?. We have plotted this 
estimate ag ains t the real pj in the spherical collapse test 
(see section 3.2) and find the maximum error to be a factor 



of ten. In practice the maximum error is dependent upon 
the problem being studied, but in general we have found the 
error to usually fall in the range of a factor of two. Note that 
the shear-free single-sided Monaghan variant would have to 
implemented as, 



Ha = 



-ap-ijCij + /3fj,jj 

Pij 



fi, 



(32) 



thereby removing the j index. Despite the estimated quan- 
tity pij, this artificial viscosity, when used in conjunction 
with the equation of motion in equation still results in a 
momentum-conserving scheme. 

In principle the V.u-based artificial viscosity should not 
suffer the problem of damping during pure shear flows, since 
the artificial viscosity only acts in compressive flows. A use- 
ful test is to supplement this artificial viscosity with the 
shear-correction term. This enables an estimate to be made 
of the extent to which the correction term under-damps due 
to SPH sampling errors. 



are given by H awley, Smarr & Wilson (1984) and Rasio & 
Shapiro (1991). In this test two regions of uniform but dif- 



ferent density gas are instantaneously brought into contact. 
If initially high density and pressure gas is to the left and the 
low density and pressure gas is to the right then a rarefac- 
tion wave will propagate left into the high density gas whilst 
a shock wave will propagate rightwards into the low den- 
sity gas. Between these there will be a contact discontinuity 
where the pressure is continuous but the density jumps. This 
behaviour is shown in Fig. |^ where the analytic solution is 
superimposed on the simulation results. 

Many authors have carried out this test in 1-dimension, 
but this is of limited use as most interesting astrophysical 
phenomena are 3-dimensional. The 1-dimensional results do 
not automatically carry over to 3-dimensions. The core of 
the SPH approach - the smoothing kernel - must be altered 
to take account of the different volume elements. Otherwise 
far too much emphasis will be placed on the central region 
in the 1-dimensional tests. Additionally, in 1-dimension, in- 
terpenetration is reduced. 

Rasio & Shapiro (1991) perform this test in 3- 
dimensions and note the large increase in numerical scatter 
and interpenetration this produces. This is partly because of 
their choice of initially randomly distributed particles which 
introduces large fluctuations in the supposedly uniform ini- 
tial density and pressure. Contrary to their assertion, this 
is not a particle distribution encountered in the course of a 
typical SPH simulation. SPH is adept at equalising density 
and pressure differences and rapidly relaxes from such an 
initial state. 



3 TEST SCENARIOS 

The SPH implementations that we examine are listed in Ta- 
ble |l| The list, whilst not exhaustive, represents a range 
of common variants for the SPH algorithm. Our motiva- 
tion is to determine the extent to which the different im- 
plementations affect behaviour in cosmological settings. In 
all cases we employ the full ft-updating scheme described in 



section 2.2 



We have chosen twelve different combinations of ar- 
tificial viscosity, symmetrization and equation of motion. 
The first is that employed by the hydra code, (Couchman, 
Pearce & Thomas 1996). We then add a shear correction 
term to this and also try a localised artificial viscosity (as 
discussed in section 2.5), resulting in three versions close 
to the TC92 formalism. We then have three versions based 



on the popular Monaghan viscosity (again see section 2.1) 
each with and without shear correction. These employ differ- 
ent ft- a vera ging schemes or kernel averaging as discussed in 
section 2.3 and are the most common schemes found in the 



literature. Finally we have three new schemes which attempt 
to combine the best features of the other implementations. 



3.1.1 Initial conditions 

We have altered our standard cubical simulation volume to 
match the geometry of a shock tube. The tube has an aspect 
ratio of 16 and all the boundaries are periodic. We take the 
high density gas to have pi — 4 and p r — 1 . The low density 
gas has p r = 1 and p r = 0.1795. In accordance with the rest 
of our tests we set y = §. 

In this test there are 4096 low density particles and 
16384 high density particles. Before starting the test itself 
each of the regions is independently evolved at constant tem- 
perature until induced random fluctuations have damped 
away. This reduces the spurious initial fluctuations men- 
tioned above and closely mimics the conditions typical of 
a cold flow in an SPH simulation. 

Note that even if the particles resided in exactly the 
correct positions to produce the theoretical curve, it is im- 
possible to reproduce the discontinuities due to the inher- 
ent smoothing in SPH estimates. The smoothed theoretical 
curves are shown in Fig. H. 



3.1.2 Shock evolution and results 



3.1 Shock tube 

The shock tube is an environment which provides an abun- 
dant source of simple tests fundament al to hydrodynami- 



cal simulations. One well-studied t est (Monaghan fc Gin- 



gold 1983; Hcrnquist fc Katz 198S; Rasio fc Shapiro 1991) 



is the Sod shock ( pocTT978 ) for which analytical solutions 



For all versions only around 35 time-steps are required 
to reach the state shown in Fig. by which time the 
shock front and the front of the rarefaction wave have 
moved around 10hi ow . The asymptotic solutions are recov- 
ered around each of the density jumps and even the small 
drop in entropy at the shock front is well modelled. Although 
the correct jump conditions are captured for the rarefaction 
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Version 


Artificial Viscosity 


symmetrization 


equation of motion 


i 


TPQ9 




i uyz 


o 
Z 


TC92-f-shcar correction 


i oyz 


i uyz 


Q 
O 


local _l uyz 


± uyz 


i uyz 


/I 


M o n ag h an 


arithmetic /ijj 


CTVTQO 

oiviyo 


r: 


Mionaghaii 


harmonic h%j 


CA/TQO 

oiviyo 


6 


Monaghan 


kernel averaging 


SM93 


7 


Monaghan+ shear correction 


arithmetic /in- 


SM93 


8 


Monaghan+ shear correction 


harmonic hij 


SM93 


9 


Monaghan+ shear correction 


kernel averaging 


SM93 


10 


Monaghan 


TC92 


TC92 


11 


Monaghan pij 


TC92 


TC92 


12 


Monaghan pij 


TC92 + kernel av 


TC92 



Table 1. Summary of the implementations examined. 'Monaghan' is the artificial viscosity of equation |26| Steinmetz-type shear correc- 
tion Steinmetz (1996) is applied where noted. The remaining terms are discussed in the text. 



wave and none of the versions introduces spurious entropy, 
the density grad ient is too shallow. Th i s effect has also been 
seen previously (Rasio & Shapiro 1991; CTP95). At the con- 



tact discontinuity the entropy jump is overshot by all the 
codes, with the problem worse for versions 1-3. The blip 
in the ideally uniform pressure at this point is due to the 
SPH s moothin g of the discontinuous density ( |Monaghan & 
Gingolr J~T983l ). The top-left panel of Fig. | shows that the 
right-moving shock front in versions 1-3 is broadened over a 
region about twice that seen for versions 4-12 which them- 
selves have about twice the width of the smoothed theoret- 
ical solution. 



3.1.3 Particle inter-penetration 

We have also used this test to study the level of flow in- 
terpenetration. Ideally there should be no interpenetration 
and the flow should remain smooth with the particles in 
well defined slabs at the end. As Fig. ^| shows for versions 1, 
2 and 3 the flow is far from smooth, with large differences 
in velocity at any one point along the shock tube. Versions 
7, 8 and 9 show much less scatter but produce noticeable 
post-shock oscillation whilst the remaining versions all pro- 
duce a very smooth flow with little velocity scatter. Due to 
the large intrinsic scatter, versions 1-3 of the code suffer 
from much greater interpenetration than the other versions. 
These effects are exacerbated in the presence of some ini- 
tial turbulence. The locally averaged nature of the viscosity 
employed by versions 1-3 allows regions to interpenetrate. 
The shear-corrected versions (7-9), despite producing just as 
sharp a density jump as the remaining versions, also perform 
poorly in the interpenetration test because interpenetrating 
streams of particles can mimic shear. 



0.4 7 

0.3 7 

0.2 7 

0.1 7 

— 

0.4 7 



0.3 

0.2 7 

0.1 7 

- 



7 Version 7 



0.4 7 -10 
3 7 Version 12 

0.2 7 
0.1 7 

: 



Figure 4. The scatter in velocity perpendicular to the shock tube 
for versions 1, 7 & 12. The x-axis is in units of the SPH search 
length in the low density gas. All the particles in the region of 
interest are shown as single points. The scatter is much higher 
and the profile is clearly different for version 1 than for either 
of the other versions. Version 12 exhibits notably less post-shock 
ringing than 7. Of the other codes, Versions 2 and 3 look like 
version 1, versions 8 and 9 look like version 7 and versions 4, 5, 
6, 10 and 11 look like version 12. 



3.1.4 Summary 

For this test the versions split into 2 main groups with the 
more accurate codes (4-12) further sub-dividing into two 
groups depending upon whether or not the viscosity is shear 
corrected. Versions 1-3 produce broader shocks and suffer 
from greater interpenetration because these codes employ 
a locally averaged method for calculating the viscosity. All 
the other versions employ a particle-particle approach to the 
viscosity calculation and hence resolve shorter scales. The 



addition of a shear-correction term to prevent transport of 
angular momentum degrades the performance of these codes 
for this particular test. Although the shock profiles are very 
similar (see Fig. |^) they produce a flow which is less smooth 
(see Fig. ^) and suffers from greater interpenetration than 
versions 4, 5, 6, 10, 11 and 12. 
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Figure 3. A comparison between the analytic result (solid abrupt line) and those from version 1 (dashed line) and version 12 (solid 
line). Also shown is the ideal result (dotted line) obtained by smoothing the theoretical result with the SPH kernel. Codes 2 and 3 match 
the profile of 1 closely, the remainder match that of 12. The x-axis is in units of the SPH smoothing length, h, in the low density gas 
and the zero point is the original interface between the high- and low-density gas. The top-left panel shows a zoom into the right-moving 
shock, the area where the largest difference between the codes is seen. 



3.2 Collapse of a spherical cloud 



3.2.1 Units and initial conditions 



This test examines the adiabatic collapse of an initially 
isothermal spherical gas cloud. We denote this test the 
'Evrard collapse' since it first appears in Evrard (1988), 
and since then has proved itself as a useful test case for 
combined SPH-gravity codes (e.g. HK89, SM93, Serna et 
al. 1996, Hultman & Kallander 1997). 



Test results are presented in normalised units. The density, 
internal energy, velocity, pressure and time are normalised 
by 3Af/47r7? 3 , GM/R, (GM/R)?, pu and (R 3 /GM)?, re- 
spectively, where R denotes the initial radius and M the 
total mass. The initial physical density distribution is given 
by, 
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30976 particles 
4776 particles 
485 particles 





Figure 5. Convergence of energy values with particle number for 
the Evrard collapse test using version 12. Energy is plotted on 
the y axis and time along the x axis in normalised units. The 
sets of curves are, from top to bottom, the thermal, kinetic, total 
and potential energies. There is comparatively little difference 
between the 4776 and 485 particle collapse because they were run 
with the same softening length. The 30796 particle run matches 
the solution calculated in SM93 very accurately. 



p{r) 



M(R) 1 
2tvR 2 r'' 



(33) 



which is achie ved by applyin g a radial stretch to an initially 
uniform grid (Evrard 19881). We prefer this configuration 



over a random one since it has sig nificantly less sampling 
erro r th an a random distribution (Huffman & Kallander 
1997j)7~|rhe internal energy of the system is chosen to be 
0.05GM / R, and the adiabatic index 7 = 5/3. The softening 
length used for each simulation is given in the caption to 
Table |. 

3.2.2 System evolution 

The evolution of the system is shown in Fig. ^. As the col- 
lapse occurs the gas is heated until the temperature of the 
core rises sufficiently to cause a 'bounce', after which a shock 
propagates outward through the gas. After the shock has 
passed through the majority of the gas, the final state of the 
system is one of virial equilibrium. Along with placing em- 
phasis on the performance of the implementations for small- 
N systems we also study convergence at larger values of N. 
A summary of the runs performed is presented in Table ^| 

3.2.3 Results from the 485 particle collapse 

With the TV considered here, the standard hydra code, and 
variants of it, do relatively poorly in this test. A lack of 
thermalization is clearly visible in Fig. |§| where the differ- 
ence in energy between versions 1-3 of the code and version 
12 are shown along the top row. These implementations all 
have thermal energy peaks that are half that of the other 
codes. The other versions perform reasonably similarly with 
minor differences being seen in the peak thermal energy and 
in the strength of post bounce oscillation (note the change 
in y-axis scale on the bottom two rows of Fig. []). 

The modified viscosity variant, 3, is slightly better at 
thermalization than versions 1 & 2 but still much worse than 



Table 2. Results of Evrard collapse test. 

Version N ste ps AE/E(xW- 3 ) AL(xl(T 7 ) 



1 


62 


6.2 


49 


485 


2 


62 


3.3 


31 


485 


3 


66 


1.1 


38 


485 


1 


91 


2.5 


11 


485 


5 


88 


1.0 


39 


485 


6 


129 


1.7 


23 


485 


7 


104 


0.8 


15 


485 


8 


95 


0.8 


65 


485 


9 


135 


1.2 


28 


485 


10 


90 


1.5 


62 


485 


11 


84 


1.5 


40 


485 


12 


103 


0.6 


29 


485 


1 


82 


6.1 


2.8 


4776 


2 


101 


6.1 


6.7 


4776 


3 


98 


1.9 


14 


4776 


1 


240 


3.1 


45 


4776 


5 


241 


3.2 


24 


4776 


6 


231 


3.5 


63 


4776 


7 


243 


3.2 


68 


4776 


8 


252 


3.1 


38 


4776 


9 


242 


3.7 


64 


4776 


10 


234 


2.9 


27 


4776 


11 


227 


2.8 


29 


4776 


12 


241 


2.8 


7.0 


4776 


1 


289 


0.3 


1.3 


30976 


3 


218 


1.9 


0.91 


30976 


1 


418 


1.4 


1.4 


30976 


6 


411 


4.5 


5.5 


30976 


11 


405 


3.6 


7.0 


30976 


12 


489 


4.7 


2.9 


30976 



Values for the 485 particle test are given at t=4.3, for the 4776 
test at t=3.4 and for the 30976 test at t=2.0. AL is measured 
in internal code units, normalising by L sta rt is not possible since 
Lstart = 0. The 485 and 4776 particle tests used a softening 
length of 0.05R, the 30976 particle tests used a softening length 
of 0.02R. 
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Figure 7. Thermal energy plot for the 485 particle collapse, com- 
paring versions 1, 3, 7 and 12. A lower peak thermalization is 
clearly visible in versions 1 and 3, whilst 7 and 12 show simi- 
lar profiles. A comparative plot of the energy difference between 
version 12 and the other codes is shown in Fig. ^ 
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Figure 6. Radial profiles for the 30976 and 485 particle collapses at t=0.8 are displayed in the left and right hand columns 
respectively. For the 30976 particle tests the density, pressure and velocity are plotted using Lagrangian bins of size 4 X 52 = 208 
particles, corresponding to 4 N grnoot f l and for the 485 particle collapses a single line connects all particles. Error bars on the 30976 
particle plots show the maximum variation within bins. The velocity plots are displaced by intervals of -0.2 vertically for easier 
interpretation. Notably versions 1 and 3 show significant differences to versions 4, 7 and 12, in particular both the density and 
pressure are underestimated (as compared to the solution calculated in SM93). Version 3 shows a better capture of the collapsing 
shock front than version 1, this effect being most evident in the pressure plot. Versions 4, 7 and 12 show excellent agreement, being 
difficult to distinguish at all radii. Version 7 shows that the inclusion of the shear-free term has little effect on the radial profile. 
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Figure 8. N=485 Evrard collapse energy difference profiles. Plots display the energy difference between version 12 and the version under consideration. Column one displays the thermal 
energy, two the kinetic energy and three the potential energy. The first row of plots compares versions 1, 2 and 3. Note the different energy scaling. The second row, displaying results 
from versions 4, 5, 6 and 7, shows the effect of the Balsara term (version 7) and also the effect of different symmetrization schemes. The third row compares version 10 and 11, thus 
indicating the effect of replacing the standard Monaghan viscosity with the single-sided version. Note that the thermal and kinetic energies in the bottom row return close to shortly 
after t=3— they do not diverge, as suggested by the graphs. Significant differences in performance, particularly for versions 1, 2 and 3, are visible for this particle number. 
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Figure 9. N=4776 Evrard collapse energy difference profiles. Plots display the difference between version 12 and the version under consideration. Column one displays the thermal 
energy, two the kinetic energy and three the potential energy. The first row of plots compares codes 1, 2 and 3. Note the different scaling of the energy. The second row, displaying results 
from codes 4, 5, 6 and 7, shows the effect of the shear-correction term. The third row demonstrates the effect of replacing the standard Monaghan viscosity (10) with the single-sided 
version (11). There is stronger convergence to a common solution than shown in the 485 particle tests. 
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any of the other versions. This indicates that the artificial 
viscosity is the primary factor in deciding the amount of ki- 
netic energy that is thermalised (as expected). The more lo- 
cal estimate of V.u used in version 3 captures the strong flow 
convergence near the bounce better than the standard esti- 
mate, leading to greater dissipation. At the other extreme 
the pairwise Monaghan viscosity, which uses the r.v trigger, 
leads to far more dissipation at bounce. It is important to 
note that the final energy values for the virialised state are 
very similar for all codes, even though their evolution is very 
different in some cases. 

Similar characteristics can be seen in the kinetic energy 
graphs. The gas in versions 1-3 develops very little kinetic 
energy. The primary cause of this is the V.w artificial vis- 
cosity which trips during the early stages of collapse (prior 
to t=0.6) and causes an increase in the thermal energy for 
all particles. This acts to decrease the compressibility of the 
gas. For all the other codes which use the Monaghan viscos- 
ity (or variant), the r.v trigger produces far less dissipation 
during the early stages of evolution, and the gas develops 
more kinetic energy. 

For the N = 485 test, Figs. ^, |^ and ^ demonstrate that 
there is no clear optimal implementation, but the general 
comparison of versions 1, 3, 7 and 12 in Fig. [?], indicates 
that some perform marginally better than others. Notable 
features of the high resolution radial PPM solution in SM93 
are a strong initial peak in the thermal energy and little post 
bounce oscillation. If we choose a model on the basis of these 
criteria then version 12 performs best, although it is difficult 
to differentiate versions 4-12 in Fig. ^. Version 12 has both 
a high initial peak and very little post bounce oscillation. It 
also conserves energy well. The shear-correction term does 
have some effect (middle row of Fig. H), in agreement with 



the observations in section 3.1. The general influence of the 



shear-correction is to increase the peak thermal energy at 
bounce and introduce slightly more post-bounce oscillation, 
although, again, this is not a significant effect. The term has 
little effect on the radial profile. 

The effect of pij replacement can be seen in the bottom 
row of Fig. |^. Versions 10 and 11 differ only by this substi- 
tution and there is very little to choose between them, the 
scatter between the kernel-averaged version 6 and version 
12 being as large. None of the implementations considered 
show poor energy conservation, and all show excellent con- 
servation of angular momentum. 



3.2.4 Effect of numerical resolution 

Increasing N to 4776 produces the expected results. The 
implementations with pairwise artificial viscosity converge 
to very similar energy profiles, see Figs. I and [u| The 10% 
spread seen in the N = 485 test is reduced to close to 1% 
and the limited scatter visible in the radial profiles is fur- 
ther reduced. The shear-correction term also has much less 
effect on the radial energy profile. For comparison, in Fig. ^ 
we show the convergence of runs performed with different 
particle number using code version 12. This plot should be 
compared with Fig. 6 of SM93. Clearly for a Monaghan- 
type viscosity the differences caused by particle number and 
softening parameter are much larger than those caused by 
the choice of SPH implementation for this range of particle 
number. 
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Figure 10. Thermal energy plot for the 4776_particle collapse, 
comparing versions 1, 3, 7 and 12. As in Fig. M, versions 1 and 
3 exhibit a lower peak thermalization although the peak value is 
marginally higher. Versions 7 and 12 again have similar profiles 
and agree very accurately on the peak thermal energy. A com- 
parative plot of the energy difference between version 12 and the 
other codes is shown in Fig. ^ 



Versions 1, 3, 4, 6, 11 and 12 were run with 30976 parti- 
cles to check for convergence of the implementations at high 
resolution. Radial profiles at t=1.4 are plotted in Fig. ^. It 
is evident from these profiles that the solutions are much 
closer than the radial profiles for the 485 particle collapse. 
However, the difference between the versions with the stan- 
dard hydra viscosity and the Monaghan variants remains 
comparatively large - a factor of two at the centre in the 
pressure and density at t=1.4. (The relatively large energy 
error is a product of our choosing a longer timestep normal- 
ization, n = 1.5 versus 1, to run the simulation in a shorter 
wall-clock time.) The profiles of the Monaghan variants all 
compare well to the radial PPM solution presented in SM93. 



3. 2. 5 Summary 

We conclude that the relatively poor shock capturing abil- 
ity of the X7.v viscosity of TC92 is a severe impediment 
to correctly calculating the evolution of this system. In 
contrast, the Monaghan viscosities (including the shear- 
corrected variants) correctly follow the evolution. All the 
Monaghan variants perform well enough in this test to be 
acceptable algorithms, and when 30976 particles are used 
in the test it is almost impossible to differentiate between 
methods (at least in the well resolved core regions). At low 
resolution (485 particles) the kernel- averaging variants pro- 
duce a slightly higher thermal peak, although some addi- 
tional ringing is visible for version 6, but not for 12. At 
medium resolution (4776 particles) convergence is stronger 
than the low resolution runs and the difference in post shock 
ringing is removed. On this basis version 12, when combined 
with its extremely fast solution time, is the preferable im- 
plementation. 



3.3 Cooling near steep density gradients 

Large density gradients occur in the gas in cosmological sim- 
ulations as a result of radiative cooling. In simulations these 
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Table 3. Cluster parameters. 



Cluster 


5 x 1O 14 M 


5 x 10 13 M o 


Rclump (kp c ) 


20 


9 


m ga3 (1O 1O M ) 


1 


0.1 


m dark (1O 1O M ) 


9 


0.9 


e (kpc) 


20 


10 



Rclum-p is the radius of the cold clump, e is the gravitational 
softening length and m ga s and m^ar*; are the mass of a gas and 
dark-matter particle respectively. 

Table 4. Cluster parameters in common. 

Cold clump Halo gas Halo dark-matter 



N 500 
T (K) 10 4 
p/p c 2 x 10 6 



4737 
2 x 10 5 - 10 8 
10 2 - 10 5 



4994 
N/A 
10 3 - 10 6 



N is the particle number, T the temperature range and p/p c the 
ratio of the density to the critical density. 

occur when cold dense knots of gas form within hot haloes 
(see section 3.5). If the smoothing radius of a hot halo par- 
ticle encompasses the cold clump then it is likely that it will 
smooth over an excess number of cold gas particles leading 
to an over-estimate of the particle's density. Consequently 
the radiative cooling for the hot particle is over-estimated, 
which allows the hot particle to cool and accrete on to the 
cold clump even though, physically, the two phases would 
be essentially decoupled. This situation should be allevi- 
ated by implementations that smooth over a fixed number 
of particles. We term this effect overcooling. It should not 
be confused with the overcooling problem (or 'cooling catas- 
trophe') in simulations of galaxy formation. 



3.3.1 Description of the halo-clump systems 

To examine this phenomenon we created core-halo systems 
each consisting of a cold clump of gas surrounded by a hot 
halo both being embedded in a dark-matter halo. The dark- 
matter and hot gas system was extracted directly from a 
cosmological simulation. The cold clump was created by ran- 
domly placing particles inside a sphere of size equal to the 
gravitational softening length and allowing this system to 
evolve to a relaxed state. The cold clump was then placed 
in the hot gas and dark matter system. Two core-halo sys- 
tems - designed to resemble galaxy clusters - were created 
to test the effect of mass and linear scale dependence, with 
total masses 5 x 10 14 M and 5 x 10 15 Mm. The parameters 
of the systems are listed in Tables |^ and Q 

Because the time-step criterion makes no reference to 
the Courant condition for cooling, it is possible that hot 
halo particles may not cool correctly as they accrete on to 
the cold clump. Tests showed that choosing the time-step 
normalization k < 0.5 was sufficient to avoid this problem. 



3.3.2 Testing the overcooling phenonomenon 

For each cluster, we prepared three experiments that exam- 
ine how the nature of the central clump changes the overall 
cooling rate. For the first experiment the central clump was 
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Figure 11. The number of cooled particles in the three experi- 
ments as a function of time. If the hot gas did not interact hydro- 
dynamically with the cold dense clump in the simulations these 
lines would overlap. 



left as a cold knot of gas, for the second it was turned into 
collisionless matter and for the third experiment the hot 
gas was allowed to become collisionless once it cooled below 
2 x 10 4 K. We denote these tests as 'standard', 'collisionless' 
and 'conversion', respectively. If the hot gas did not inter- 
act with cold gas then these tests would all give the same 
result. The number of particles cooled over time for these 
three tests for the 5 x 10 14 M cluster is shown in Fig. 

The behaviour at early times (< 2 x 10 9 yrs) is dom- 
inated by a sudden rise in the number of cold particles. 
This is due to the hot gas within 2h of the dense clump 
responding to its sudden introduction. In the standard test, 
gas particles comprise the dense clump and hence the den- 
sities calculated for the hot gas rise suddenly, causing rapid 
cooling. For both the collisionless and conversion tests, the 
dense clump is collisionless; the hot-gas densities rise only in 
response to contraction of the halo about the clump. Over- 
cooling does not become significant in the collisionless test 
until a sufficient number of gas particles (approximately 40) 
have cooled and contracted in the central region to provide 
a 'seed' clump. Once formed, this seed clump permits the 
rapid overcooling of the rest of the hot gas particles in its 
immediate vicinity. Consequently, a comparable number of 
particles are cooled during this period as in the standard 
test. By construction the conversion test never forms this 
seed clump and hence overcooling is never initiated. Cool- 
ing occurs only by the contraction of the hot gas halo. 

At later times (> 2 x 10 9 yrs), the greater number of 
gas particles in the central region for the standard test (600 
versus 200 for the collisionless test) creates a larger gas- 
density gradient and hence more efficient overcooling. The 
accretion is fed by a contracting h alo. It is this quasi-steady 
state which is examined in Section 3.3.3. For the conversion 



test, the increasing central mass density, coupled with the 
absence of a central gas clump to provide pressure support, 
leads to a progressively increasing cooling rate which is not 
directly related to the overcooling phenomenon. 

The a-physical drop in temperature of the 'hot' gas par- 
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Figure 13. The amount of overcooled gas for two clusters of dif- 
ferent mass. The y-axis is the percentage excess of cooled gas for a 
'standard' run relative to the amount cooled in the corresponding 
'collisionless' run. 



tides when the cold clump falls within twice their smoothing 
length is clearly illustrated by comparing the temperature 
profile of the gas produced in the standard test with that 
produced by the conversion test (Figs. |l| (a) and (b)). The 
conversion test produces a profile that approximately re- 
sembles that to be expected in the absence of hot and cold 
gas phase interaction. For this test the gas is approximately 
isothermal except in the core, where the increased density 
has caused the gas to cool - this is clearly different to the 
standard test. Note that in the conversion test the central 
core is more extended than in the standard test, but remains 
within a sphere smaller than the hot gas smoothing radius. 

The behaviour of the gas at the interface between the 
gas phases is illustrated in Fig. [l^. For the standard test the 
smoothing process forces the density to rise very abruptly 
from the halo to the core, whilst for the conversion test the 
lack of a cold gas core removes this imperative. Consequently 
in the standard test, particles outside the dense core, but 
within 2hhot, have a high density and cooling rate. Thus once 
particles fall within 2hhot an abrupt temperature decrease 
results as the cooling rate increases - as seen in panel (a). In 
the conversion test the flatter density profile does not lead 
to excessively high cooling rates, and there is no abrupt 
temperature drop. 

The cooling rate of the gas is also affected by the virial 
temperature of the halo gas. The overcooling for the two 
different clusters are compared in Fig. uA 



cance of the observed variations we note that four different 
realisations of the same test produced a significant variation 
of up to 40 cooled particles after 2 x 10 9 yrs. We expect, how- 
ever, that for a given realisation the general trend amongst 
the SPH variants would be the same. 

The primary source of variation in the overcooling rate 
among the implementations of SPH is the viscosity (Fig. |l4| 
a). The V.v-based artificial viscosities produce cooling rates 
that are essentially indistinguishable, while the Monaghan 
variants lead to a significantly greater cooling (about 50% 
more). This difference is probably related to the V.u variants 
providing somewhat greater pressure support in the core (see 
section 3.2). The inclusion of a shear-correction term in the 
artificial viscosity (Fig. [w| b) has little effect on the cooling 
rate - as expected. 

Since the overcooling effect is caused by the large dif- 
ference in kernel sizes associated with the hot halo particles 
and the cold clump particles, it might be expected that the 
symmetrization has a role to play in determining the cooling 
rate. Consider the ^-averaging schemes: the arithmetic mean 
is limited to having a minimum value of hi arge /2, while the 
harmonic mean is zero if any particle interacts with another 
particle having h — 0. In practice, as is shown in Figs. |l4| 
c) and d), there is little difference between all symmetriza- 
tion schemes, with the exception of version 11. This version 
combines a pure gather kernel, with the single-sided Mon- 
aghan artificial viscosity. This result is surprising in view 
of the comparatively 'normal' results for versions 10, which 
differs in terms of the artificial viscosity, and 12 which differ 
in terms of the symmetrization. 



3.3.4 Summary 

All the versions of SPH we have tested exhibit overcool- 
ing and this effect should be seen as generic to the method 
itself. SPH will always experience difficulties modelling ar- 
bitrarily steep density gradients. The only implementation 
that stands out as performing poorly is 11 which couples 
a one-sided implementation of Monaghan artificial viscosity 
with the TC92 symmetrization procedure. When the TC92 
symmetrization is supplemented with kernel averaging the 
problem is removed. 



3.4 Drag 

There is concern that the over-merging problem encountered 
in N-body simulations of clusters of galaxies is exacerbated 
in simulations which use SPH (see Frenk et al. 1996). Exces- 
sive drag on small knots of gas within a hot halo will cause 
the knots to spiral inward into regions of stronger tidal forces 
where they may be disrupted (e.g., Moore et al. 1996). In 
this section we model a cold dense clump moving through 
a hot halo and investigate if the problem is sensitive to the 
particular SPH implementation employed. 



3.3.3 Results of the different SPH implementations 

Experiments were run with all 12 different SPH implemen- 
tations. The collisionless tests produced cooling rates which 
were essentially identical. We report on implementation sen- 
sitivities for the standard test. In order to assess the signifi- 



3.4-1 Drag test model systems 

To cover a variety of infall speeds we examine the decelera- 
tion of a knot of cold gas in three velocity regimes: Mach 2, 
Mach 1, and Mach 1/3. The Mach 2 and Mach 1 tests differ 
in terms of the speed of the cold knot ('fast' versus 'slow') 
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Figure 12. Temperature and density profiles for the standard test (a, c) and conversion test (b, d). The mean smoothing radii for both 
the hot and cold particles are represented by the vertical lines. The line for h shows the point internal to which the bulk of the kernel 
weighting is accumulated. The 2h line gives the outer limit of the smoothing. The inner pair of lines (left) are for the cold gas, while 
the outer pair are for the hot gas. In the lower panel, the cold particles are collisionless, and hence have no smoothing length. For these 
data, the initial conditions of the 5 X 10 14 Mq mass cluster were used. In the conversion test the stars are plotted at the temperature to 
which they had cooled just before conversion to star particles. 



and not the temperature of the hot gas. The Mach 1/3 test 
uses the same clump velocity as the Mach 1 test, but is per- 
formed in hotter ('very hot') gas. Table JE] gives the details 
of the cold clump and hot gas phases. Clump characteristics 
are selected to loosely emulate a poorly resolved galaxy with 
no dark matter, while the hot gas media are typical of the 
intracluster medium. 

The hot gas was prepared from an initially random 
placement of particles, and then allowed to relax to a sta- 
ble state. The cold clump was created by randomly placing 
particles within a sphere of radius equal to the gravitational 
softening length. The cold clump was allowed to relax in 
the same manner as the hot gas, before combining the two 
systems. 

The Jeans length, Rj, for the hot gas phases is suffi- 
ciently large to ensure stability even in the presence of the 
perturbation from the cold clump. Consequently, dynamical 
friction should not be important. This conclusion was con- 



firmed by passing a collisionless cold clump through the hot 
medium - it experienced negligible deceleration. 

The box length, 5 Mpc, was chosen so that the cold 
clump was well separated from its images (arising from the 
periodic boundary conditions employed) and would move 
across the bo x on ly once without encountering its own wake. 
As in Section 3.3, an appropriate value of the time-step nor- 



malization parameter, k, was found. For these tests, a value 
of k — 1.0 is used. 



3.4-2 Expected deceleration 

An expected rate of deceleration may be approximated by 
considering a disc sweeping through a hot medium, collect- 
ing all matter it encounters. This would represent a maxi- 
mum expected rate of deceleration, ignoring dynamical fric- 
tion. The solution for the velocity, V, of such a system is 
given by V(t) = l/(t — ti), where I is a characteristic length 
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Figure 14. Variation with time in the number of cold particles (T < 10 5 K) due to: a) different implementations of the artificial 
viscosity. The viscosity given by TC92, both in original forms (dashed line) and in the more localized implementation (dot-dashed line) 
lead to a similar amount of cooling as the artificial viscosity given by Monaghan (dotted line). However, the one-sided version of the 
Monaghan viscosity (upper solid line) does produce an extra excess of cold particles. For comparison, the amount of cooling produced 
in the run with an initially collisionless dense core is given by the solid line, b) the presence of a Balsara term in the artificial viscosity, 
c) the symmetrization scheme. Those implementations using arithmetically averaged values of the smoothing lengths, h, produce the 
mean cooling curve given by the dashed line. The mean curved produced by those that use a harmonically averaged value of h give the 
dot-dashed line. If the kernels themselves are averaged, the dotted line is the result, d) the symmetrization of TC92 and its variant. 
The mean of the rates given in panel c) (dashed line) is just slightly lower than the rates produced by the code version which uses the 
symmetrization procedure of TC92 as well as that of its more localized variant, both of which are essentially the same. These data are 
for the 5 X 10 14 M n cluster. 



given by I = M/2nR 2 p g and ti is a characteristic time-scale 
given by ti = t — l/V - Here, M is the mass of the disc at 
the start, R is the radius of the disc, and p g is the density 
of the gas through which the disc is travelling. The clump 
starts with velocity V at time t . For the tests that use the 
slow clump, this estimate implies the final velocity should be 
400 km/s. For the fast clump, the final clump velocity should 
be 670 km/s. These crude estimates indicate that hydrody- 
namical forces should indeed be important for the parame- 
ters being considered. 



3.4-3 Results of the SPH variants 

Four separate realizations of the same initial conditions were 
evolved with the same version of the test code to look for 
variation due to randomness in the initial conditions. There 
is variation on the order of 10 km/s between the runs for 
the Mach 2 and Mach 1/3 scenarios, 4 km/s for the Mach 1 
scenario. 

The cold clump size varies between implementations. It 
is 2-3 times larger for the V .v viscosity variants. However, 
since the knot size is still much less than the smoothing 
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Table 5. The characteristics of the cold clumps and the hot media 
used in the drag tests. 





Slow cold clump 


Fast cold clump 


P/Pc 


1000 


1000 


T (K) 


10 4 


10 4 


R (kpc) 


50 


50 


N 


100 


100 


m (10 9 M Q ) 


1.7 


1.7 


Vo (km/s) 


500 


1000 




Hot gas 


Very hot gas 


P/Pc 


10 


10 


T (K) 


10 7 


10 8 


TV 


13000 


13000 


m (10 9 M Q ) 


1.7 


1.7 


V s (km/s) 


500 


1500 


Rj (Mpc) 


6 


18 



Given are the overdensity, p/p c (/iioo = l)i the temperature, T, 
the radius of the cold clump, R, the number of particles in the 
medium, N, the mass resolution of the medium, m, the initial 
velocity of the cold clump, V , the speed of sound in the hot 
medium, V s , and the Jeans length for the hot medium, Rj. The 
simulation volume in all cases is (5 Mpc) 3 . The 'fast cold clump' 
was used in the Mach 2 runs in combination with the 'hot gas'. 
The Mach 1 runs used the 'slow cold clump' embedded in the 'hot 
gas'. The Mach 1/3 runs used the 'slow cold clump' in the 'very 
hot gas'. 



Table 6. The relative final velocities of the cold clumps. 



Version 


Mach 2 


Mach 1 


Mach 1/3 


1 


0.999 ± 0.002 


1.004 ± 0.009 


1.24 ± 0.10 


2 


0.991 ± 0.006 


1.003 ± 0.012 


1.29 ± 0.08 


3 


0.914 ±0.006 


1.009 ± 0.008 


1.20 ± 0.16 


1 


1.018 ± 0.005 


1.009 ± 0.008 


0.78 ± 0.07 


5 


1.025 ± 0.009 


1.004 ± 0.009 


0.84 ± 0.07 


6 


0.984 ± 0.008 


0.926 ± 0.013 


0.67 ± 0.05 


7 


1.064 ± 0.007 


1.060 ± 0.004 


1.26 ± 0.02 


8 


1.045 ± 0.007 


1.062 ± 0.008 


1.34 ±0.01 


9 


1.028 ± 0.006 


0.979 ± 0.006 


1.12 ± 0.05 


10 


0.957 ±0.005 


0.951 ± 0.007 


0.58 ± 0.07 


11 


0.956 ± 0.004 


0.955 ± 0.006 


0.62 ± 0.06 


12 


1.018 ± 0.003 


1.038 ± 0.009 


1.06 ± 0.10 



Given is the mean relative velocity of the cold clumps over the 
final 0.5 X 10 9 yrs normalized by the mean velocity of all the cold 
clumps in that velocity regime. 

radius for these particles (at least a factor of three), the total 
clump size, after smoothing, is approximately the same in 
all cases. 

Compared with V.u viscosity, Monaghan viscosity in 
both the symmetric and single-sided forms leads to an in- 
crease in the damping of the velocity of the clump when 
used with the TC92 symmetrization (Table ^ and Fig. |l5|) . 
However, when Monaghan viscosity is used with TC92 sym- 
metrization, supplemented by kernel averaging, the deceler- 
ation becomes comparable to the V.u versions. The more 
localized estimate of the V.w viscosity does little except in 
the Mach 2 set of runs, for which it increases the drag to 
match that of the Monaghan viscosity. The inclusion of a 
shear-correction term reduces the drag in the Mach 1/3 case 
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Figure 15. Variation of the cold-clump velocity with artificial 
viscosity type (no shear-correction term is included). In each 
panel the different lines distinguish different viscosity implemen- 
tations: standard HYDRA viscosity (TC92), solid; localized TC92 
viscosity, dashed; standard Monaghan viscosity, dash-dot; single- 
sided Monaghan viscosity, dotted. The panels are, from top to 
bottom; Mach 2, Mach 1 and Mach 1/3. 

as well as in the Mach 1 case when the clump velocity has 
dropped below Mach 0.8. 

Use of either the arithmetic or harmonic average for hij 
produces less drag then any other symmetrization method 
(Fig. [l6]) except the version with TC92 symmetrization com- 
bined with kernel averaging. On their own, kernel averaging 
and the TC92 symmetrization produce marginally higher 
deceleration, particularly at supersonic speeds. 

3.4-4 Summary 

The tests favour (but cannot distinguish between) the har- 
monic and arithmetic averages. The shear-correction term 
lowers the drag at sub-sonic speeds. The Monaghan viscos- 
ity coupled with the TC92 symmetrization performs poorly. 

3.5 Cosmological simulation 

In this test we simulate a common astrophysical problem: 
the formation of knots of cold, dense gas within a cosmolog- 
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Figure 16. Variation of the cold-clump velocity with h- 
symmctrization. The lines are: arithmetic average, solid; har- 
monic, dashed; kernel averaging, dot-dash; TC92 symmetrization 
(version 10), dotted. The panels are, from top to bottom; Mach 
2, Mach 1 and Mach 1/3. 



ical volume. For this problem we are interested in recovering 
accurate positions and masses for the objects. However, the 
resolution is such that no internal information (such as spiral 
structure or radial density profiles) can be recovered. 



3.5.1 Initial Conditions 

The simulations presented here were of an firj = 1, standard 
cold dark matter universe with a box size of 10/i -1 Mpc. We 
take h = 0.5 throughout this section, equivalent to a Hubble 
constant of SOkms" 1 Mpc" 1 . The baryon fraction, Qi, was 



set from nucleosynthesis constraints, Qth 2 = 0.015 (Copi ei 
1995) and we assume a constant gas metallicity of Q.5Zq. 



Identical initial conditions were used in all cases, allowing a 
direct comparison to be made between the objects formed. 

The initial fluctuation amplitude was set by requir- 
ing that the model produce the same number-density of 
rich clusters as observed today. To achieve this we take 
erg = 0.6, the present-day linear rms fluctuation on a scale 
of 8ft" 1 Mpc (Eke, Cole & Frenk 1996, Vianna & Liddle 
1996). Each model began with 32 3 dark matter particles 
each of mass 1.58 x 10 10 Mq and 32 3 gas particles each of 



Table 7. Properties 


of the cosmolog 


ical test 


runs by 


version. 




Version 


AT 

1 v steps 


Hours 


AT 


M ■ 

big 


ga 


hot 


f 

J cold 


i 


6060 


34.7 


30 


1.01 


0.17 


0.37 


0.45 


2 


6348 


35.0 


34 


1.07 


0.19 


0.36 


0.45 


3 


6532 


43.7 


41 


1.08 


0.20 


0.34 


0.46 


1 


6984 


65.1 


49 


2.41 


0.30 


0.37 


0.33 


5 


6113 


60.5 


48 


2.34 


0.27 


0.40 


0.33 


6 


6769 


57.4 


56 


Z.Oo 


U.o4 


U.o4 


U.oo 


7 


7562 


63.4 


47 


2.43 


0.29 


0.34 


0.36 


8 


6782 


93.6 


47 


2.38 


0.27 


0.36 


0.36 


9 


7688 


109.7 


19 


2.40 


0.32 


0.37 


0.31 


10 


6858 


60.5 


50 


2.69 


0.33 


0.33 


0.34 


11 


6522 


40.4 


48 


2.43 


0.31 


0.35 


0.34 


12 


6205 


38.5 


51 


2.28 


0.30 


0.35 


0.34 



Listed are the number of steps taken to reach 2 = 0, the number 
of hours required on a Sun Ultra II 300 workstation, the number 
of groups of more than 50 cold particles found at the endpoint, 
the mass of the largest clump (in units of 10 12 Mq), the fraction 
of the gas in galaxies, the fraction of gas above 10 5 K and the 
fraction of the gas that remains diffuse and cold (all at z = 0). 



mass 1.01 x 1O 9 M0, smaller than the critical mass derived 
by Steinmetz & White (1997) required to prevent 2-body 
heating of the gaseous component by the heavier dark mat- 
ter particles. The simulations were started at redshift 19. We 
employ a comoving Plummer softening of 10/i - kpc, which 
is typical for modern cosmological simulations but still larger 
than required to accurately simulate the dynamics of galax- 
ies in dense environments. This test case is identical to that 
extensively studied by Kay et al. (1998) who used it to exam- 
ine the effect of changing numerical and physical parameters 
for a fixed SPH implementation. 



3.5.2 Extraction of glob properties 

The gas is effectively in three disjoint phases. There is a 
cold, diffuse phase which occupies the dark-matter voids and 
therefore most of the volume. A hot phase occupies the dark- 
matter halos and at this resolution is typically above 10 5 K. 
Finally, there is a cold, dense phase consisting of tight knots 
of gas typically at densities several thousand times the mean 
and at temperatures close to 10 4 K. The relative proportions 
of the gas in each of these phase s is g iven in Table |?]. We 
follow Evrard, Summers & Davis ( 1994 ) in defining a cooled 



knot of particles as a 'glob' because the resolution is such 
that they can hardly be termed a galaxy. The properties 
of the globs are calculated by first extracting all the parti- 
cles which are simultaneously below a temperature of 10° K 
and at densities above 180 times the mean and then run- 
ning a friend- of- friends group-finder with a maximum link- 
ing length, b, of 0.08 times the mean interparticle separation 
of the dark matter. In practice the object set obtained is in- 
sensitive to the choice of b because the globs are typically 
disjoint, tightly bound clumps. The cumulative multiplic- 
ity function for the different implementations is shown in 
Fig. 0. 
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Figure 17. The multiplicity functions for different versions. The 
six implementations plotted span the range of outcomes as, in all 
cases, the addition of a shear-correction to the viscosity makes 
little difference. Versions 10 and 11 produce very similar results 
to version 12, whilst 1, 2 and 3 produce smaller objects than the 
other versions. 
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Figure 18. Comparison of object masses in different implemen- 
tations. The masses of objects found by the group-finder in ver- 
sion 12 of the code are compared to the masses found for the 
corrsponding objects in other versions. Versions 1-3 of the code 
all produce objects of about half the mass and many smaller ob- 
jects are missing (because in these runs they fall below the reso- 
lution limit of around 50 particles). As for Fig. jl?] the addition of 
a shear-correction makes little difference and versions 10 and 11 
of the code produce very similar masses to version 12. 



3.5.3 Results of cosmological test 

In all cases the lar gest object has 'overcooled' in the sense de- 



scribed in section 3.3, It is much too massive to be expected 
in a simulation of this size and is only present because gas 
within the hot halo has its cooling rate enhanced by the very 
high-density gas contained in the globs. 

A distinct difference can be seen in the morphologies 
of small objects formed by versions 1-3 compared to those 
formed by 4-12. Versions 1-3 produce spherical objects since 



the X7.v viscosity used does not damp random orbital motion 
within the softening radius. Versions 4-12 produce disc-like 
objects as a result of the effective dissipation provided by 
the pairwise trigger and conservation of angular momentum. 
Both spherical and disc objects are of size approximately 
equal to the softening length. We do not expect merging 
to play a significant role in this simulation due to the low 
particle-number in the majority of globs. 

The major discriminant between the versio ns i s the dif- 
ferent artificial viscosities. As shown in section 3.1, versions 
1, 2 and 3 produce broader shock fronts because the viscos- 
ity employed is a locally averaged quantity, whereas for all 
the other versions the viscosity is calculated on a pairwise 
basis. The pairwise viscosity shock-heats the gas more effi- 
ciently and leads to a larger fraction of hot gas in versions 
4-12, (see Table 0). 

The fraction of matter present in globs and the number 
of groups with N > 50 is clearly lower for versions 1-3 than 
for other versions. The lower mass-fraction in globs is a com- 
bination of both the smaller number of groups found above 
the threshold and versions l-3_producing lighter objects. As 



was demonstrated in section 3.2, the V.u viscosity produces 
a shallow collapse, with much less dissipation. When col- 
lapse occurs in an object that has approximately N smoot h 
particles, there will be virtually no shock heating and the 
gas particles will free-stream within the shallow potential 
well. In simulations with cooling the V.w viscosity provides 
a marginally higher pres sure support than the Monaghan 
viscosity (see section |3.6[ ) which can be sufficient to prevent 
collapse of surrounding material. At low resolution this re- 
sults in an object not achieving N > 50, whilst at higher 
resolution the object has lower mass. 

Of the Monaghan variants 4-6, version 5 has the lowest 
fraction of matter in the glob phase and the highest fraction 
of hot gas. Fig. jl8] shows that it also produces systemati- 
cally lighter objects. Version 6 has the highest fraction of 
matter in the glob phase, the lowest fraction of hot gas and 
tends to produce the heaviest objects. These results are due 
to the symmetrization scheme causing the artificial viscos- 
ity to produce different amounts of dissipation. Similarly, 
for versions 10-12 the fraction of hot gas can be traced to 
the amount of d issip ation. These results are consistent with 
those in section 3.6, where they are discussed in detail. The 
trend for Monaghan- type viscosities is distinct: versions that 
produce more dissipation form lighter objects as the hot halo 
gas is heated to higher temperatures where the cooling time 
is longer. 

In section 3.1 it was shown that the shear-correction 
term is less able to capture shocks and consequently pro- 
duces lower shock heating. For the /i-averaging implemen- 
tations (4, 5), the hot gas fraction is reduced upon adding 
shear correction, which agrees with the shock tube result. 
For the kernel- averaged version this is not the case - the 
hot gas fraction inc reases. This result is probably not signif- 
icant; in section |3~6| versions 4-6 all show reduced dissipation 
upon including the shear-correction term. 



3.5.4 Summary 

Any of the versions discussed in this paper could be used ef- 
fectively for this problem. Differences in the amount of gas 
in each of the hot, cold and glob phases are produced, which 
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Table 8. Results of rotating cloud collapse test. 



Ver 


N s tep 


AE/E(xlO~ 3 ) 


AL/L(xlO~ 4 ) 


Qpeak 


fpeak 


1 


1461 


2.S 


+0.23 


0.001 


0.095 


2 


1627 


1.7 


+0.98 


0.001 


0.118 


3 


1588 


1.7 


+3.8 


0.024 


0.212 


4 


1806 


2 8 


+9.9 


0.023 


0.159 


5 


1703 


1.8 


+9.2 


0.049 


0.215 


6 


1834 


1.8 


+10. 


0.009 


0.124 


7 


2087 


4.4 


-6.2 


0.015 


0.128 


8 


2002 


4.4 


-6.9 


0.030 


0.164 


9 


2074 


4.8 


-5.7 


0.004 


0.095 


10 


1789 


2.2 


+10. 


0.002 


0.074 


11 


1748 


3.1 


+7.4 


0.023 


0.143 


12 


1705 


2.9 


+6.6 


0.036 


0.175 



With the exception of the peak variables values are given at 
t=256. AL/L, the fractional change in angular momentum, is 
measured to be positive for a loss of momentum and is expected 
to be zero. Q pea k is the maximum value of the thermal energy 
(as a fraction of the initial total mechanical energy) and f pea k is 
the peak fraction of gas shocked to high temperature. 

can be explained in terms of the amount of dissipation pro- 
duced by each scheme. The V.u-based viscosities produce 
objects of spherical morphology, while the Monaghan vari- 
ants produce objects with disc morphology. 
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Figure 20. Comparison of morphology for N = 1736 run under 
varying h m i n , at t = 128. Each panel is 100 kpc across, and grav- 
itational softening was set at 2 kpc. Implementation 12 was used 
to run the simulations. The right-hand panel clearly shows more 
structure on scales near and below the gravitational softening 
length. 

momentum in the initial conditions a 'ring' of dark matter 
is thrown off. Swing amplification causes transitory spiral 
features early in the evolution which are later replaced by 
spiral structure that persists for a number of rotations. If 
shocked gas is developed during the collapse it forms a halo 
around the disc. 



3.6 Rotating cloud collapse 

A standard test problem for galaxy - formation codes is pre- 
sented in Navarro & White ( 1995 ). In this test a cloud of 
dark matter and gas is set in solid-body rotation. Gravita- 
tional collapse combined with radiative cooling leads to a 
cool, centrifugally supported gaseous disc. 



3.6.1 Initial conditions 

The initial radius of the gas cloud is chosen to be 100 kpc, 
and the total mass (dark matter and gas) is 10 12 Mq. The 
spin parameter, A, is set to, 



A = 



|£||g| 1/2 
GM 5 / 2 



0.1, 



(34) 



where L is the angular momentum, E the binding energy, 
and M the mass. The baryon fraction is set to 10%, and 
the cooling function for primordial- abund ance gas is inter- 
polated from Sutherland & Dopita (1993). 

For our tests we consider simulations with N — 2 x 1736 
particles. The gravitational softening length is set at 2 kpc 
for both dark matter and gas particles. This is different from 
previous authors who have set the dark-matter softening to 
be 5 kpc and the gas softening to 2 kpc. As a result we form 
a smaller central dark-matter core. Times are quoted in the 



units of Navarro & White ( 1993) (4.7 x 10 fa yr). One rotation 



period at the half-mass radius corresponds to approximately 
170 timesteps. 

The resulting evolution of the system is shown in 
Fig. |l9[ and test results summarized in Table [| Radiative 
cooling during the collapse causes the gas to form a fiat 
disc. The dark matter virialises quickly after collapse, leav- 
ing a tight core. Because of the large amount of angular 



3.6.2 Non-implementation-specific results 

We have found marginally different results for our codes 
when compared with other work. This is due to two factors. 
Firstly, using a 2 kpc softening for the dark-matter particles 
has a significant effect on the final morphology. The 2-body 
interaction between gas and dark matter is much stronger 
than would be expected if the dark matter had a longer 
softening length. Secondly, most of the particles in the disc 
have an h value close to h m i n which in turn sets a signif- 
icant limit on the minimum mass of a clump that may be 
resolved. We have run a simulation with h m in = 0.05e to see 
the effect of this. Fig. ^o| shows a comparison of the simula- 
tion run with the smaller h m in to the standard hmin = 0.5e 
simulation. Far more structure is evident on scales close to 
the gravitational softening length, which must be viewed as 
being unphysical since at this scale the gravitational forces 
are severely softened. 

Since the circular velocity is calculated from [GM(< 
r )/ r ] 1 ^ 2 an d the dark matter has the dominant mass contri- 
bution we expect little difference among the rotation curves 
for the different implementations. In Fig. ^ we plot the ro- 
tation curves for four different implementations. Apart from 
a visibly lower central mass concentration for version 1 there 
is comparatively little difference. 



3.6.3 Implementation- specific results 

Before the disc has formed (prior to t=128) versions 1-3 
have an extended gas halo compared with the remaining 
versions. The halo for version 1 is as much as 40% larger 
than those for versions 4-12. In Fig. ^ we compare the gas 
structure of version 1 to that of version 12. The source of the 
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Figure 19. Evolution of gas and dark matter in 2 X 1736 particle collapse. The results for version 10 are plotted, which produces 
little shocked gas during collapse. The morphological evolution of the system agrees well with previous work, with minor differences 
being attributable to differing initial conditions. The results presented here preserve symmetry above and below the equatorial 
plane for longer than seen in other work, this may be a consequence of the excellent momentum conservation exhibited by grid 
based gravity solvers and our smoothing over all particles within 2h m i n . 
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Figure 21. Comparison of gas morphology for 2 X 1736 particle collapse. Time is t=256 and axis scales are in kpc. Clearly different 
implementations exhibit different spiral structures, indicating that at this resolution the structures are poorly defined. 
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Figure 22. Rotation curves for four different implementations. 
The shear correction in version 9 produces a higher rotational ve- 
locity at half the disc radius. This is because the outward trans- 
port of angular momentum is reduced, thereby reducing the in- 
ward movement of the half mass radius. 
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Figure 23. Thermal energy, divided by the total initial mechani- 
cal energy, for 5 different implementations. A marginal change in 
the artificial viscosity can produce a 25-fold change in the peak 
total thermal energy. This is indicative of the test sitting at the 
edge of the Rees-Ostriker (1977) criterion. 



extended halo is the V.w artificial viscosity, which acts to in- 
crease the local pressure. This is seen in the early rise in the 
thermal energy for version 1 in Fig. ^3| The more local es- 
timate used in version 3 produces less pressure support and 
a smaller halo. The Monaghan viscosity does not provide 
pressure support as the pairwise r.v term is very small. The 
different artificial viscosities also lead to different disc mor- 
phologies. The V.w viscosity fails to damp collapse along the 
z-axis sufficiently and allows far more interpenetration than 
the Monaghan viscosity leading to thicker discs in versions 
1-3. 

The angular momentum losses in Table ^ show a notice- 
able trend. For most codes AL/L is small and positive (by 
definition indicating a loss of angular momentum). However 
the shear-corrected Monaghan variants show an increase in 
the angular momentum of the system. However, since the 
magnitude of the angular momentum is approximately the 
same as that of the other codes, we do not place strong signif- 
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Figure 24. Fraction of gas mass above 30,000 K. This gives an 
approximate measure of the amount of gas which is shocked. The 
amount of shocked gas is extremely sensitive to the artificial vis- 
cosity implementation. Compare versions 10 and 11 which differ 
only by the pij substitution, the peak amount of shocked gas is 
different by a factor of two. 



icance on this result. Version 2 also has the shear-correction 
term, but we attribute the similar performance to version 1 
as being due to the low amount of dissipation produced by 
the V.w viscosity. 

Examining the thermal energy during collapse yields 
very interesting results. Fig. shows a plot of the thermal 
energy of the cloud versus time. As a fraction of the total 
energy the thermal energy makes a small contribution be- 
cause the baryon fraction is only 10%. However, the relative 
differences in thermal energy between versions can be sig- 
nificant. This situation is analogous to the differences seen 
in t he kinetic energy in the Evrard collapse test (see section 
3.2). For this test it is important to note that the differences 
in the thermal energy arise from the amount of shocked gas 
present in the simulation, which is determined by the artifi- 
cial viscosity. This is demonstrated in the comparison plot of 
the fraction of gas above 30000 K, shown in Fig. |24|. A small 
change in the artificial viscosity can have a very significant 
change in the amount of shocked gas. If one considers the 
effect of changing the ft-averaging scheme from the arith- 
metic mean to the harmonic mean, then the artificial vis- 
cosity will be higher in most situations. This is because the 



term used to prevent divergences, 0.01ft;.,- in the denomina- 
tor of equation is now always smaller, and hence can lead 
to larger values of TLij. This explains why version 5 has so 
much shocked gas. Similarly, the pij replacement in versions 
11 and 12 leads to a larger as the (hi/hj) 3 pi replacement 
systematically tends to underestimate the value of pj, and 
hence we get more shocked gas. We note, however, that the 
effect is only noticeable in cases of extreme density contrast, 
e.g., halo particles just above a cold gaseous disc. 

It is also evident that a change in the symmetrization 
procedure can have a significant effect, codes 4 and 6 have 
differing amounts of shocked gas. This fact suggests that 
the halo gas in this collapse problem mus t sit at the edge of 
the Rees-Ostriker ( Rees fc Ostriker 1977 ) cooling criterion, 
namely that the free fall time is approximately equal to the 
gas cooling time. We have checked this by running simula- 
tions with masses a factor of five higher and two lower. The 
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leads to a disc that is unstable to perturbations ( |Navarro~& 
White 1993). The limit at which this o ccurs is set by the 
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Figure 25. Comparison of gas halo size at t = 128 for versions 
1 and 12. Axis scales are given in kpc. The gas halo for version 1 
is clearly larger. 

lower mass system produces less (12% by mass, compared 
to 18%) shocked gas, whilst the high mass system produces 
a very large (75%) amount of shocked gas. In view of these 
results and that the thermal energy is a very small fraction 
of the total, we do not place strong emphasis on the dif- 
ferences in shocked gas betwe en th e versions. Comparison 
of the res ults o f Serna et al. (1996) and those of Navarro 
& White (1993) confirms this conclusion as the amount of 
shocked gas in the former differs visibly from the latter. 



3.6.4 Summary 

Setting aside the differences in the amount of shocked gas 
among implementations, there is comparatively little varia- 
tion among the final results. There are differences in mor- 
phology: versions 1, 2 and 3 have a thicker disc and, during 
the initial collapse, show a more extended gas halo. Both 
effects can be traced to the V.i> viscosity. The potential and 
kinetic energies during collapse, however, are all similar, due 
to the dominance of the dark matter. 



3.7 Disc stability 

Disc stability is a critica l issu e in galaxy formation. It is now 
widely known (Balsara (1995), Navarro & Steinmetz (1997), 
for example) that the standard Monaghan viscosity intro- 
duces spurious angular momentum transport (as opposed 
to the physical transfer of angular momentum that occurs 
in differentially rotating discs). This spurious transfer can 
have a significant effect of the development of a small N 
object, as the angular momentum may be tra nsported to 
the outer edge of the disc in a few rotations (Navarro & 
Steinmetz 1997| ). 

In this section we compare the growth of the half- 
angular-momentum radius (half- AM radius) of the disc, the 
half-mass radius and ratio of the two to determine which 
SPH implementation is least susceptible to this problem. 
A similar investig ation was first performed by (Navarre 
& Steinmetz 1997). For initial conditions we use the disc 



formed by version 9 during the rotating cloud collapse (see 



section 



3.6). We cut out the central 100 kpc diameter region 



Toomre stability criterion (Toomre 1964), Q, where 



Q 



(35) 



3.36GE ' 

oy is the velocity dispersion, v the epicycle frequency and E 
the surface density. 

3. 7. 1 General evolution properties 

Given the comparatively quiet disc environment the mor- 
phological evolution of the different versions were compara- 
tively similar. The only noticeable difference could be seen 
in the disc thickness, which for the V.u viscosity variants (1, 
2 and 3), was much larger than that of the rest of the ver- 
sions. This effect was also observed in section 3.6, although 



and evolved this for a sufficient time for relaxation of the sys- 
tem as well as for transport of angular momentum. We did 
not run a higher resolution test as increasing the resolution 



here the gas has 'diffused' away from an initially thin disc. 



3.7.2 Implementation-specific results 

The evolution of the half-AM radii are plotted in Fig. |26[ 
The figure indicates that the codes fall into two groups, with 
one group suffering a stronger decay of the half-mass radius 
and an associated growth of the half-AM radius. The other 
group, which shows less decay, consists of versions 1-3, 7-9. 
The artificial viscosities for this group are the V.u variants 
(1-3) and the shear-corrected Monaghan version (7-9). For 
this group the half-AM radius does not change significantly 
during the simulation which corresponds to approximately 
30 rotations and 5000 time-steps. The half-mass radius de- 
cays by approximately 50%, leading to a similar reduction 
in the half-mass to half-AM radii ratio. 

Within this first group of codes there is a sub-division 
determined by the disc thickness at the end of the simula- 
tion. The V.t> variants all have a thick disc, due to the failure 
of the algorithm to adequately damp convergent motions in 
the z-direction. For the shear-corrected Monaghan variants 
this is not a concern. 

The second group, comprising versions 4-6, 10-12 (all 
Monaghan viscosity variants), shows significant outward 
transport of the angular momentum, and by the end of the 
simulation the half radius has increased by approximately 
50%. There is also a larger decay in the half-mass radius, it 
being 50% greater than the decay seen for the other group 
of codes. These two results contribute to make the half-mass 
to half-AM radii ratio decay to only 25% of its initial value 
at the beginning of the simulation. This result confirms that 
seen in Navarro & Steinmetz (1997) where it was shown that 
the shear correction significantly reduces angular momen- 
tum transport in disc simulations. Our results show a faster 
increase in transport, but this is probably due to the sim- 
ulations being dissimilar; Navarro & Steinmetz place a disc 
in a predefined halo, and then evolve the combined system. 



3.7.3 Summary 

Whilst preserving the half-mass to half-AM radii ratio as 
well as the shear-corrected scheme, the V.u viscosity is not a 
significant improvement. The large increase in disc thickness 
and loss of definition more than outweighs the improvement 
in the decay of the half-mass radius. Both the Monaghan and 



© 1998 RAS, MNRAS 000, 



28 R. J. Thacker et al. 



- 

'-5 



02 

a 

X 

B 
S 
a 



- 

a 

< 



D2 

< 
ac 

o 




0.05 



200 300 400 500 600 700 800 900 1000 

Time 



Figure 26. The half-mass radius, half-AM radius and the ratio of the two. Versions 1 and 9 are best, producing a negligible 
increase in the angular momentum half-mass radius and a comparatively slow reduction in the half-mass radius. The Monaghan 
variants that are not shear-corrected perform worst, transporting the angular momentum rapidly, resulting in a fast decrease of the 
half- mass radius. Radii are plotted in code units (1 unit equals 800 kpc). 
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shear-corrected artificial viscosities maintain the disc struc- 
ture. On the basis of these tests any of the shear-corrected 
implementations (7-9) is to be preferred. 



4 CONCLUSIONS 

We have presented a series of tests designed to determine 
the differences in performance of various SPH implemen- 
tations in scenarios common in simulations of cosmological 
hierarchical clustering. Special attention was paid to how 
the codes perform for small- N problems. A summary of our 
findings is presented in Table bl 

Our principle conclusions follow. 

(i) We recommend schemes that use the Monaghan vis- 
cosity supplemented with the shear-correction. Of those 
methods that do not use the shear-correction, version 12 
is preferable because of its high speed and accuracy (and 
see vii below). 

(ii) Several implementations introduce programming dif- 
ficulties, such as changing the neighbour search from a 
gather process to a hybrid gather-scatter. This is especially 
problematic in the adaptive refinements where different sym- 
metrization schemes lead to different choices for which par- 
ticles to place in a refinement. The TC92 symmetrization is 
by far the easiest to program in this respect. 

(iii) The choice of artificial viscosity is the primary factor 
in determining code performance. The equation of motion 
and particle symmetrization schemes have only a secondary 
role, albeit significant for small- TV. In particular, the artifi- 
cial viscosity used in the current implementation of hydra, 
and variants of it, produce a large amount of scatter in local 
variables, such as the velocity field and temperature, and 
also lead to less thermalisation. These characteristics indi- 
cate that the relative performance of an artificial viscosity 
is determined by its effective resolution. Viscosities which 
use an estimate of V.u to determine whether the viscosity 
should be applied will always be less able to capture strong 
flow convergence and shocks than those which use the pair- 
wise r.v trigger. 

(iv) Instabilities inherent in simple smoothing-length up- 
date algorithms can be removed. By using a 'neighbour- 
counting kernel' and a weighted update average, stability 
can be increased without requiring an expensive exact cal- 
culation of the correct h to yield a constant number of neigh- 
bours. 

(v) We strongly agree with the conclusion presented in 
SM93 that to accurately calculate local physical variables 
in dynamically evolving systems, at least 10 4 particles are 
required. The difference in morphologies observed in the ro- 
tating cloud collapse clearly indicates that the belief SPH 
can accurately predict galaxy morphologies with as few as 
1000 particles is overly optimistic. This implies that in stud- 
ies of galaxy formation, where the internal dynamics define 
morphology, an object may only be considered well resolved 
if it contains a minimum of 10 4 particles. 

(vi) The introduction of a shear-corrected viscosity leads 
to reduced shock capturing, although the effect is small and 
primarily visible in the local velocity field. We have also con- 
firmed the results of other authors, namely that the shear- 
corrected viscosity does indeed reduce viscous transport of 
angular momentum. 



(vii) It is possible to implement a scheme (11 and 12) 
which uses the Monaghan artificial viscosity, that does not 
require the precomputation of all density values before solv- 
ing the hydrodynamic equation of motion. Further, the re- 
sulting implementation conserves momentum and energy to 
an extremely high accuracy. This implementation removes 
the need to compute the list of neighbours twice or alterna- 
tively the need to store it for every particle. Hence it is sig- 
nificantly faster than other schemes. Additionally, although 
not done in this work, the shear-correction can be incorpo- 
rated with little extra effort. 
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Table 9. Qualitative summary of the strengths and weaknesses of each implementation. We categorize each version using a y/ to indicate 
preferable performance and an X to indicate inferior performance. + and - signs differentiate between similarly grouped implementations. 
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algorithm. In high-resolution simulations (2 X 128 3 particles, for example) versions 4—9 may well be significantly slower. Section numbers 
from which these conclusions are drawn are shown at the top of the columns. 
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